home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / zgeev.f < prev    next >
Text File  |  1996-07-19  |  13KB  |  384 lines
  1.       SUBROUTINE ZGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR,
  2.      $                  WORK, LWORK, RWORK, INFO )
  3. *
  4. *  -- LAPACK driver routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       CHARACTER          JOBVL, JOBVR
  11.       INTEGER            INFO, LDA, LDVL, LDVR, LWORK, N
  12. *     ..
  13. *     .. Array Arguments ..
  14.       DOUBLE PRECISION   RWORK( * )
  15.       COMPLEX*16         A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
  16.      $                   W( * ), WORK( * )
  17. *     ..
  18. *
  19. *  Purpose
  20. *  =======
  21. *
  22. *  ZGEEV computes for an N-by-N complex nonsymmetric matrix A, the
  23. *  eigenvalues and, optionally, the left and/or right eigenvectors.
  24. *
  25. *  The right eigenvector v(j) of A satisfies
  26. *                   A * v(j) = lambda(j) * v(j)
  27. *  where lambda(j) is its eigenvalue.
  28. *  The left eigenvector u(j) of A satisfies
  29. *                u(j)**H * A = lambda(j) * u(j)**H
  30. *  where u(j)**H denotes the conjugate transpose of u(j).
  31. *
  32. *  The computed eigenvectors are normalized to have Euclidean norm
  33. *  equal to 1 and largest component real.
  34. *
  35. *  Arguments
  36. *  =========
  37. *
  38. *  JOBVL   (input) CHARACTER*1
  39. *          = 'N': left eigenvectors of A are not computed;
  40. *          = 'V': left eigenvectors of are computed.
  41. *
  42. *  JOBVR   (input) CHARACTER*1
  43. *          = 'N': right eigenvectors of A are not computed;
  44. *          = 'V': right eigenvectors of A are computed.
  45. *
  46. *  N       (input) INTEGER
  47. *          The order of the matrix A. N >= 0.
  48. *
  49. *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
  50. *          On entry, the N-by-N matrix A.
  51. *          On exit, A has been overwritten.
  52. *
  53. *  LDA     (input) INTEGER
  54. *          The leading dimension of the array A.  LDA >= max(1,N).
  55. *
  56. *  W       (output) COMPLEX*16 array, dimension (N)
  57. *          W contains the computed eigenvalues.
  58. *
  59. *  VL      (output) COMPLEX*16 array, dimension (LDVL,N)
  60. *          If JOBVL = 'V', the left eigenvectors u(j) are stored one
  61. *          after another in the columns of VL, in the same order
  62. *          as their eigenvalues.
  63. *          If JOBVL = 'N', VL is not referenced.
  64. *          u(j) = VL(:,j), the j-th column of VL.
  65. *
  66. *  LDVL    (input) INTEGER
  67. *          The leading dimension of the array VL.  LDVL >= 1; if
  68. *          JOBVL = 'V', LDVL >= N.
  69. *
  70. *  VR      (output) COMPLEX*16 array, dimension (LDVR,N)
  71. *          If JOBVR = 'V', the right eigenvectors v(j) are stored one
  72. *          after another in the columns of VR, in the same order
  73. *          as their eigenvalues.
  74. *          If JOBVR = 'N', VR is not referenced.
  75. *          v(j) = VR(:,j), the j-th column of VR.
  76. *
  77. *  LDVR    (input) INTEGER
  78. *          The leading dimension of the array VR.  LDVR >= 1; if
  79. *          JOBVR = 'V', LDVR >= N.
  80. *
  81. *  WORK    (workspace/output) COMPLEX*16 array, dimension (LWORK)
  82. *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  83. *
  84. *  LWORK   (input) INTEGER
  85. *          The dimension of the array WORK.  LWORK >= max(1,2*N).
  86. *          For good performance, LWORK must generally be larger.
  87. *
  88. *  RWORK   (workspace) DOUBLE PRECISION array, dimension (2*N)
  89. *
  90. *  INFO    (output) INTEGER
  91. *          = 0:  successful exit
  92. *          < 0:  if INFO = -i, the i-th argument had an illegal value.
  93. *          > 0:  if INFO = i, the QR algorithm failed to compute all the
  94. *                eigenvalues, and no eigenvectors have been computed;
  95. *                elements and i+1:N of W contain eigenvalues which have
  96. *                converged.
  97. *
  98. *  =====================================================================
  99. *
  100. *     .. Parameters ..
  101.       DOUBLE PRECISION   ZERO, ONE
  102.       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  103. *     ..
  104. *     .. Local Scalars ..
  105.       LOGICAL            SCALEA, WANTVL, WANTVR
  106.       CHARACTER          SIDE
  107.       INTEGER            HSWORK, I, IBAL, IERR, IHI, ILO, IRWORK, ITAU,
  108.      $                   IWRK, K, MAXB, MAXWRK, MINWRK, NOUT
  109.       DOUBLE PRECISION   ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
  110.       COMPLEX*16         TMP
  111. *     ..
  112. *     .. Local Arrays ..
  113.       LOGICAL            SELECT( 1 )
  114.       DOUBLE PRECISION   DUM( 1 )
  115. *     ..
  116. *     .. External Subroutines ..
  117.       EXTERNAL           DLABAD, XERBLA, ZDSCAL, ZGEBAK, ZGEBAL, ZGEHRD,
  118.      $                   ZHSEQR, ZLACPY, ZLASCL, ZSCAL, ZTREVC, ZUNGHR
  119. *     ..
  120. *     .. External Functions ..
  121.       LOGICAL            LSAME
  122.       INTEGER            IDAMAX, ILAENV
  123.       DOUBLE PRECISION   DLAMCH, DZNRM2, ZLANGE
  124.       EXTERNAL           LSAME, IDAMAX, ILAENV, DLAMCH, DZNRM2, ZLANGE
  125. *     ..
  126. *     .. Intrinsic Functions ..
  127.       INTRINSIC          DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN, SQRT
  128. *     ..
  129. *     .. Executable Statements ..
  130. *
  131. *     Test the input arguments
  132. *
  133.       INFO = 0
  134.       WANTVL = LSAME( JOBVL, 'V' )
  135.       WANTVR = LSAME( JOBVR, 'V' )
  136.       IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
  137.          INFO = -1
  138.       ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
  139.          INFO = -2
  140.       ELSE IF( N.LT.0 ) THEN
  141.          INFO = -3
  142.       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  143.          INFO = -5
  144.       ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
  145.          INFO = -8
  146.       ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
  147.          INFO = -10
  148.       END IF
  149. *
  150. *     Compute workspace
  151. *      (Note: Comments in the code beginning "Workspace:" describe the
  152. *       minimal amount of workspace needed at that point in the code,
  153. *       as well as the preferred amount for good performance.
  154. *       CWorkspace refers to complex workspace, and RWorkspace to real
  155. *       workspace. NB refers to the optimal block size for the
  156. *       immediately following subroutine, as returned by ILAENV.
  157. *       HSWORK refers to the workspace preferred by ZHSEQR, as
  158. *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
  159. *       the worst case.)
  160. *
  161.       MINWRK = 1
  162.       IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
  163.          MAXWRK = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 )
  164.          IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN
  165.             MINWRK = MAX( 1, 2*N )
  166.             MAXB = MAX( ILAENV( 8, 'ZHSEQR', 'EN', N, 1, N, -1 ), 2 )
  167.             K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'ZHSEQR', 'EN', N, 1,
  168.      $          N, -1 ) ) )
  169.             HSWORK = MAX( K*( K+2 ), 2*N )
  170.             MAXWRK = MAX( MAXWRK, HSWORK )
  171.          ELSE
  172.             MINWRK = MAX( 1, 2*N )
  173.             MAXWRK = MAX( MAXWRK, N+( N-1 )*
  174.      $               ILAENV( 1, 'ZUNGHR', ' ', N, 1, N, -1 ) )
  175.             MAXB = MAX( ILAENV( 8, 'ZHSEQR', 'SV', N, 1, N, -1 ), 2 )
  176.             K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'ZHSEQR', 'SV', N, 1,
  177.      $          N, -1 ) ) )
  178.             HSWORK = MAX( K*( K+2 ), 2*N )
  179.             MAXWRK = MAX( MAXWRK, HSWORK, 2*N )
  180.          END IF
  181.          WORK( 1 ) = MAXWRK
  182.       END IF
  183.       IF( LWORK.LT.MINWRK ) THEN
  184.          INFO = -12
  185.       END IF
  186.       IF( INFO.NE.0 ) THEN
  187.          CALL XERBLA( 'ZGEEV ', -INFO )
  188.          RETURN
  189.       END IF
  190. *
  191. *     Quick return if possible
  192. *
  193.       IF( N.EQ.0 )
  194.      $   RETURN
  195. *
  196. *     Get machine constants
  197. *
  198.       EPS = DLAMCH( 'P' )
  199.       SMLNUM = DLAMCH( 'S' )
  200.       BIGNUM = ONE / SMLNUM
  201.       CALL DLABAD( SMLNUM, BIGNUM )
  202.       SMLNUM = SQRT( SMLNUM ) / EPS
  203.       BIGNUM = ONE / SMLNUM
  204. *
  205. *     Scale A if max element outside range [SMLNUM,BIGNUM]
  206. *
  207.       ANRM = ZLANGE( 'M', N, N, A, LDA, DUM )
  208.       SCALEA = .FALSE.
  209.       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
  210.          SCALEA = .TRUE.
  211.          CSCALE = SMLNUM
  212.       ELSE IF( ANRM.GT.BIGNUM ) THEN
  213.          SCALEA = .TRUE.
  214.          CSCALE = BIGNUM
  215.       END IF
  216.       IF( SCALEA )
  217.      $   CALL ZLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
  218. *
  219. *     Balance the matrix
  220. *     (CWorkspace: none)
  221. *     (RWorkspace: need N)
  222. *
  223.       IBAL = 1
  224.       CALL ZGEBAL( 'B', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
  225. *
  226. *     Reduce to upper Hessenberg form
  227. *     (CWorkspace: need 2*N, prefer N+N*NB)
  228. *     (RWorkspace: none)
  229. *
  230.       ITAU = 1
  231.       IWRK = ITAU + N
  232.       CALL ZGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
  233.      $             LWORK-IWRK+1, IERR )
  234. *
  235.       IF( WANTVL ) THEN
  236. *
  237. *        Want left eigenvectors
  238. *        Copy Householder vectors to VL
  239. *
  240.          SIDE = 'L'
  241.          CALL ZLACPY( 'L', N, N, A, LDA, VL, LDVL )
  242. *
  243. *        Generate unitary matrix in VL
  244. *        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
  245. *        (RWorkspace: none)
  246. *
  247.          CALL ZUNGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
  248.      $                LWORK-IWRK+1, IERR )
  249. *
  250. *        Perform QR iteration, accumulating Schur vectors in VL
  251. *        (CWorkspace: need 1, prefer HSWORK (see comments) )
  252. *        (RWorkspace: none)
  253. *
  254.          IWRK = ITAU
  255.          CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VL, LDVL,
  256.      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
  257. *
  258.          IF( WANTVR ) THEN
  259. *
  260. *           Want left and right eigenvectors
  261. *           Copy Schur vectors to VR
  262. *
  263.             SIDE = 'B'
  264.             CALL ZLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
  265.          END IF
  266. *
  267.       ELSE IF( WANTVR ) THEN
  268. *
  269. *        Want right eigenvectors
  270. *        Copy Householder vectors to VR
  271. *
  272.          SIDE = 'R'
  273.          CALL ZLACPY( 'L', N, N, A, LDA, VR, LDVR )
  274. *
  275. *        Generate unitary matrix in VR
  276. *        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
  277. *        (RWorkspace: none)
  278. *
  279.          CALL ZUNGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
  280.      $                LWORK-IWRK+1, IERR )
  281. *
  282. *        Perform QR iteration, accumulating Schur vectors in VR
  283. *        (CWorkspace: need 1, prefer HSWORK (see comments) )
  284. *        (RWorkspace: none)
  285. *
  286.          IWRK = ITAU
  287.          CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VR, LDVR,
  288.      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
  289. *
  290.       ELSE
  291. *
  292. *        Compute eigenvalues only
  293. *        (CWorkspace: need 1, prefer HSWORK (see comments) )
  294. *        (RWorkspace: none)
  295. *
  296.          IWRK = ITAU
  297.          CALL ZHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, W, VR, LDVR,
  298.      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
  299.       END IF
  300. *
  301. *     If INFO > 0 from ZHSEQR, then quit
  302. *
  303.       IF( INFO.GT.0 )
  304.      $   GO TO 50
  305. *
  306.       IF( WANTVL .OR. WANTVR ) THEN
  307. *
  308. *        Compute left and/or right eigenvectors
  309. *        (CWorkspace: need 2*N)
  310. *        (RWorkspace: need 2*N)
  311. *
  312.          IRWORK = IBAL + N
  313.          CALL ZTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
  314.      $                N, NOUT, WORK( IWRK ), RWORK( IRWORK ), IERR )
  315.       END IF
  316. *
  317.       IF( WANTVL ) THEN
  318. *
  319. *        Undo balancing of left eigenvectors
  320. *        (CWorkspace: none)
  321. *        (RWorkspace: need N)
  322. *
  323.          CALL ZGEBAK( 'B', 'L', N, ILO, IHI, RWORK( IBAL ), N, VL, LDVL,
  324.      $                IERR )
  325. *
  326. *        Normalize left eigenvectors and make largest component real
  327. *
  328.          DO 20 I = 1, N
  329.             SCL = ONE / DZNRM2( N, VL( 1, I ), 1 )
  330.             CALL ZDSCAL( N, SCL, VL( 1, I ), 1 )
  331.             DO 10 K = 1, N
  332.                RWORK( IRWORK+K-1 ) = DBLE( VL( K, I ) )**2 +
  333.      $                               DIMAG( VL( K, I ) )**2
  334.    10       CONTINUE
  335.             K = IDAMAX( N, RWORK( IRWORK ), 1 )
  336.             TMP = DCONJG( VL( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
  337.             CALL ZSCAL( N, TMP, VL( 1, I ), 1 )
  338.             VL( K, I ) = DCMPLX( DBLE( VL( K, I ) ), ZERO )
  339.    20    CONTINUE
  340.       END IF
  341. *
  342.       IF( WANTVR ) THEN
  343. *
  344. *        Undo balancing of right eigenvectors
  345. *        (CWorkspace: none)
  346. *        (RWorkspace: need N)
  347. *
  348.          CALL ZGEBAK( 'B', 'R', N, ILO, IHI, RWORK( IBAL ), N, VR, LDVR,
  349.      $                IERR )
  350. *
  351. *        Normalize right eigenvectors and make largest component real
  352. *
  353.          DO 40 I = 1, N
  354.             SCL = ONE / DZNRM2( N, VR( 1, I ), 1 )
  355.             CALL ZDSCAL( N, SCL, VR( 1, I ), 1 )
  356.             DO 30 K = 1, N
  357.                RWORK( IRWORK+K-1 ) = DBLE( VR( K, I ) )**2 +
  358.      $                               DIMAG( VR( K, I ) )**2
  359.    30       CONTINUE
  360.             K = IDAMAX( N, RWORK( IRWORK ), 1 )
  361.             TMP = DCONJG( VR( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
  362.             CALL ZSCAL( N, TMP, VR( 1, I ), 1 )
  363.             VR( K, I ) = DCMPLX( DBLE( VR( K, I ) ), ZERO )
  364.    40    CONTINUE
  365.       END IF
  366. *
  367. *     Undo scaling if necessary
  368. *
  369.    50 CONTINUE
  370.       IF( SCALEA ) THEN
  371.          CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, W( INFO+1 ),
  372.      $                MAX( N-INFO, 1 ), IERR )
  373.          IF( INFO.GT.0 ) THEN
  374.             CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, W, N, IERR )
  375.          END IF
  376.       END IF
  377. *
  378.       WORK( 1 ) = MAXWRK
  379.       RETURN
  380. *
  381. *     End of ZGEEV
  382. *
  383.       END
  384.